############################################################
# Replication: ALL Figures (1–7)
# Output dir:  ~/Dropbox/Issue voting Chile/09_replication/output
############################################################

rm(list = ls())

# ---- Packages ----
library(haven)
library(foreign)
library(readstata13)
library(dplyr)
library(forcats)
library(ggplot2)
library(ggeasy)
library(cregg)
library(gridExtra)
library(Hmisc)
library(lfe)
library(broom)
library(rdrobust)
library(rdd)

# ---- Paths ----
base_dir <- "~/Dropbox/Issue voting Chile"
outdir   <- file.path(base_dir, "09_replication/output")
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)

# ---- Helper: always save with white background ----
save_fig <- function(filename, plot, width = 6, height = 4.5, dpi = 300) {
  ggsave(filename, plot = plot, width = width, height = height, dpi = dpi, bg = "white")
}

# ---- Shared palette ----
dark_palette <- c("#990000", "#003399", "#006600", "black")

############################################################
# Figure 1: Ideological Identification over time (CEP)
############################################################
cep <- haven::read_dta(file.path(base_dir, "09_replication/data/cep_all.dta"))

left     <- aggregate(left ~ encuesta_a, data = cep, FUN = mean)
right    <- aggregate(right ~ encuesta_a, data = cep, FUN = mean)
center   <- aggregate(center ~ encuesta_a, data = cep, FUN = mean)
none_pol <- aggregate(none_pol ~ encuesta_a, data = cep, FUN = mean)

left$ideology     <- "Left"
right$ideology    <- "Right"
center$ideology   <- "Center"
none_pol$ideology <- "No ideology"

left     <- dplyr::rename(left,     Percentage = left)
right    <- dplyr::rename(right,    Percentage = right)
center   <- dplyr::rename(center,   Percentage = center)
none_pol <- dplyr::rename(none_pol, Percentage = none_pol)

ideology_all <- rbind(left, right, center, none_pol)
ideology_all$ideology   <- as.factor(ideology_all$ideology)
ideology_all$encuesta_a <- as.factor(ideology_all$encuesta_a)

ideology_plot <- ggplot(
  ideology_all,
  aes(x = encuesta_a, y = Percentage, color = factor(ideology), group = factor(ideology))
) +
  geom_line() +
  labs(x = "Years", y = "Ideological Identification (%)") +
  scale_color_manual(values = dark_palette, name = NULL) +
  scale_x_discrete(c(1994, 2005, 2007, 2009, 2011, 2013, 2015, 2017, 2019)) +
  theme_minimal() +
  theme(
    axis.text.x  = element_text(size = 12),
    axis.text.y  = element_text(size = 12),
    axis.title.y = element_text(size = 14),
    legend.text  = element_text(size = 14)
  )

save_fig(file.path(outdir, "figure1_ideology.png"), ideology_plot, width = 7.5, height = 4.8)

############################################################
# Figure 2: Partisan Identification over time (CEP)
############################################################
party    <- aggregate(party ~ encuesta_a, data = cep, FUN = mean)
no_party <- aggregate(no_party ~ encuesta_a, data = cep, FUN = mean)

party$partisanship    <- "Any Party"
no_party$partisanship <- "No Party"

party    <- dplyr::rename(party,    Percentage = party)
no_party <- dplyr::rename(no_party, Percentage = no_party)

party_all <- rbind(party, no_party)
party_all$partisanship <- as.factor(party_all$partisanship)

party_plot <- ggplot(
  party_all,
  aes(x = encuesta_a, y = Percentage, color = factor(partisanship), group = factor(partisanship))
) +
  geom_line() +
  labs(x = "Years", y = "Partisan Identification (%)") +
  scale_color_manual(values = dark_palette, name = NULL) +
  scale_x_continuous(breaks = seq(1994, 2023, 2)) +
  theme_minimal() +
  theme(
    axis.text.x  = element_text(size = 12),
    axis.text.y  = element_text(size = 12),
    axis.title.y = element_text(size = 14),
    legend.text  = element_text(size = 14)
  )

save_fig(file.path(outdir, "figure2_partisanship.png"), party_plot, width = 7.5, height = 4.8)

############################################################
# Figures 3 & 4: Panel stacked bars (2021 vs 2023)
############################################################
panel_long <- read_dta(file.path(base_dir, "09_replication/data/panel_long.dta"))

# Keep panel_all==1 & waves 2 or 3
panel_all <- subset(panel_long, panel_all == 1 & (wave == 2 | wave == 3))

# Optional: normalize ideology label
panel_all$ideology <- forcats::fct_recode(panel_all$ideology, "No Ideology" = "Independent")

## Figure 3: Ideology composition by wave
pct_ideol <- panel_all %>%
  group_by(wave, ideology) %>%
  summarise(count = n(), .groups = "drop_last") %>%
  mutate(percentage = count / sum(count)) %>%
  ungroup() %>%
  filter(!is.na(ideology))

fig3_ideology_wave <- ggplot(pct_ideol, aes(x = wave, y = percentage, fill = ideology)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = paste0(round(percentage * 100, 1), "%")),
            position = position_stack(vjust = 0.5),
            size = 4, color = "white", show.legend = FALSE) +
  labs(x = "Year", y = "") +
  scale_fill_manual(values = dark_palette, name = NULL) +
  scale_x_continuous(breaks = c(2, 3), labels = c("2021", "2023")) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2)) +
  theme_minimal() +
  theme(
    axis.text.x   = element_text(size = 14, color = "black"),
    axis.title.x  = element_text(size = 16),
    axis.text.y   = element_blank(),
    legend.text   = element_text(size = 12),
    panel.border  = element_rect(color = "gray", fill = NA, size = 1)
  )
save_fig(file.path(outdir, "figure3_ideology_wave.png"), fig3_ideology_wave, width = 6.2, height = 5)

## Figure 4 (left): Immigration by wave
pct_mig <- panel_all %>%
  group_by(wave, immigration) %>%
  summarise(count = n(), .groups = "drop_last") %>%
  mutate(percentage = count / sum(count)) %>%
  ungroup()

fig4_left_mig <- ggplot(pct_mig, aes(x = wave, y = percentage, fill = immigration)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = paste0(round(percentage * 100, 1), "%")),
            position = position_stack(vjust = 0.5),
            size = 4, color = "white", show.legend = FALSE) +
  labs(x = "Year", y = "") +
  scale_fill_manual(values = c("#990000", "#006600", "#003399", "black"), name = NULL) +
  scale_x_continuous(breaks = c(2, 3), labels = c("2021", "2023")) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2)) +
  theme_minimal() +
  theme(
    axis.text.x   = element_text(size = 14, color = "black"),
    axis.title.x  = element_text(size = 16),
    axis.text.y   = element_blank(),
    legend.text   = element_text(size = 12),
    panel.border  = element_rect(color = "gray", fill = NA, size = 1)
  )
save_fig(file.path(outdir, "figure4_left_immigration_wave.png"), fig4_left_mig, width = 6.2, height = 5)

## Figure 4 (right): Feminism by wave
pct_fem <- panel_all %>%
  group_by(wave, A47_wave) %>%
  summarise(count = n(), .groups = "drop_last") %>%
  mutate(percentage = count / sum(count)) %>%
  ungroup()

pct_fem <- pct_fem %>%
  mutate(A47_wave = forcats::fct_recode(A47_wave,
                                        "Yes" = "Sí",
                                        "Don't Know" = "No sé"))

fig4_right_fem <- ggplot(pct_fem, aes(x = wave, y = percentage, fill = A47_wave)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = paste0(round(percentage * 100, 1), "%")),
            position = position_stack(vjust = 0.5),
            size = 4, color = "white", show.legend = FALSE) +
  labs(x = "Year", y = "") +
  scale_fill_manual(values = dark_palette, name = NULL) +
  scale_x_continuous(breaks = c(2, 3), labels = c("2021", "2023")) +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2)) +
  theme_minimal() +
  theme(
    axis.text.x   = element_text(size = 14, color = "black"),
    axis.title.x  = element_text(size = 16),
    axis.text.y   = element_blank(),
    legend.text   = element_text(size = 12),
    panel.border  = element_rect(color = "gray", fill = NA, size = 1)
  )
save_fig(file.path(outdir, "figure4_right_feminism_wave.png"), fig4_right_fem, width = 6.2, height = 5)

# =========================
# Conjoint: Figures 5 & 6
# =========================

############################################################
# Figures 5 & 6 – Conjoint (mm)
# Output: ~/Dropbox/Issue voting Chile/09_replication/output
############################################################

rm(list = ls())

# ---- Packages ----
library(haven)
library(dplyr)
library(forcats)
library(ggplot2)
library(ggeasy)
library(cregg)

# ---- Paths ----
base_dir <- "~/Dropbox/Issue voting Chile"
outdir   <- file.path(base_dir, "09_replication/output")
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)

# Always save with a white background (avoid grey when transparent)
save_fig <- function(filename, plot, width = 5, height = 5.4, dpi = 300) {
  ggsave(filename, plot = plot, width = width, height = height, dpi = dpi, bg = "white")
}

# ---- Load data (prefer replication path, fallback to clean_data) ----
w2_path_main <- file.path(base_dir, "09_replication/data/wave2_conjoint.dta")
w2_path_alt  <- file.path(base_dir, "01_data/clean_data/wave2_conjoint.dta")
wave2_conjoint <- if (file.exists(w2_path_main)) read_dta(w2_path_main) else read_dta(w2_path_alt)

# ---- Prepare attributes as factors and harmonize labels ----
atr <- c("a_1","a_2","a_3","a_4","a_5","a_6")
wave2_conjoint[, atr] <- lapply(wave2_conjoint[, atr], factor)

wave2_conjoint <- wave2_conjoint %>%
  mutate(
    Ideology    = a_1,
    Gender      = a_2,
    Age         = a_3,
    Feminism    = a_4,
    Immigration = a_5,
    Crime       = a_6
  )

# Clean levels consistently
wave2_conjoint$Ideology <- fct_collapse(
  wave2_conjoint$Ideology,
  Left  = c("Left"),
  Right = c("Right")
)
wave2_conjoint$Ideology <- factor(wave2_conjoint$Ideology, levels = c("Left","Right"))

wave2_conjoint$Feminism <- fct_collapse(
  wave2_conjoint$Feminism,
  `Non-Feminist` = c("Non-Feminist"),
  Feminist       = c("Feminist")
)
wave2_conjoint$Feminism <- factor(wave2_conjoint$Feminism, levels = c("Non-Feminist","Feminist"))

wave2_conjoint$Crime <- fct_recode(
  wave2_conjoint$Crime,
  "Less Punitive" = "No Harsher Punishment",
  "More Punitive" = "Harsher Punishment"
)
wave2_conjoint$Crime <- factor(wave2_conjoint$Crime, levels = c("Less Punitive","More Punitive"))

wave2_conjoint$Immigration <- fct_recode(
  wave2_conjoint$Immigration,
  "New Restrictions" = "New Immigration Restrictions"
)
wave2_conjoint$Immigration <- factor(
  wave2_conjoint$Immigration,
  levels = c("No Restrictions","New Restrictions")
)

# Construct subgroup flags if missing
if (!("left_pro_mig" %in% names(wave2_conjoint))) {
  wave2_conjoint$left_pro_mig <- as.integer(
    (wave2_conjoint$left == 1) &
      (wave2_conjoint$a6e %in% c("Ni de acuerdo ni en desacuerdo","En desacuerdo","Muy en desacuerdo"))
  )
}
if (!("right_anti_mig" %in% names(wave2_conjoint))) {
  wave2_conjoint$right_anti_mig <- as.integer(
    (wave2_conjoint$right == 1) &
      (wave2_conjoint$a6e %in% c("De acuerdo","Muy de acuerdo"))
  )
}

dl     <- function(x) droplevels(x)
fix_ci <- function(df) {
  if (!("conf.low" %in% names(df)) || !("conf.high" %in% names(df))) {
    stopifnot(all(c("estimate", "std.error") %in% names(df)))
    df <- df %>%
      mutate(conf.low  = estimate - 1.96 * std.error,
             conf.high = estimate + 1.96 * std.error)
  }
  df
}

# Generic plotter: black, thinner points/CIs (size = 0.55)
mm_draw <- function(mm_df, title, order_levels = NULL) {
  df <- fix_ci(mm_df)
  if (!is.null(order_levels)) {
    df <- df %>% dplyr::filter(level %in% order_levels)
    df$level <- factor(df$level, levels = order_levels) |> droplevels()
  }
  
  ggplot(df, aes(x = estimate, y = level)) +
    geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
                    color = "black", size = 0.12) +
    geom_vline(xintercept = 0.3) +
    scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2)) +
    labs(title = title, x = NULL, y = NULL) +
    theme_bw() + guides(color = "none") +
    easy_center_title() +
    theme(
      plot.title  = element_text(size = 16, color = "black"),
      axis.text.x = element_text(size = 14, color = "black"),
      axis.text.y = element_text(size = 12, color = "black")
    )
}

# Extract Ideology & Immigration rows from an mm() result
select_ideo_mig <- function(mm_obj) {
  mm_obj %>%
    filter(feature %in% c("Ideology", "Immigration")) %>%
    mutate(level = as.character(level))
}

# -------------------------
# Figure 5 – NON‑INTERACTED
# -------------------------

# Left & Pro‑Immigration (non‑interacted)
left_pro <- dl(subset(wave2_conjoint, left_pro_mig == 1))
mm_left_all <- mm(
  left_pro,
  choice_clean2 ~ Ideology + Gender + Age + Feminism + Immigration + Crime,
  id = ~ numericalid
)
mm5_left <- select_ideo_mig(mm_left_all)

fig5_left <- mm_draw(
  mm5_left,
  title = "Left Pro-Immigration",
  order_levels = c("Left","Right","No Restrictions","New Restrictions")
)
save_fig(file.path(outdir, "mm_left_mig_non.png"), fig5_left, width = 5, height = 5.4)

# Right & Anti‑Immigration (non‑interacted)
right_anti <- dl(subset(wave2_conjoint, right_anti_mig == 1))
mm_right_all <- mm(
  right_anti,
  choice_clean2 ~ Ideology + Gender + Age + Feminism + Immigration + Crime,
  id = ~ numericalid
)
mm5_right <- select_ideo_mig(mm_right_all)

fig5_right <- mm_draw(
  mm5_right,
  title = "Right Anti-Immigration",
  order_levels = c("Left","Right","No Restrictions","New Restrictions")
)
save_fig(file.path(outdir, "mm_right_mig_non.png"), fig5_right, width = 5, height = 5.4)

# ----------------------------------------------
# Figure 6 – Ideology × Immigration (interacted)
# ----------------------------------------------

mm_interaction <- function(dat, title) {
  dat <- dl(dat)
  dat$profile1 <- interaction(dat$Ideology, dat$Immigration, sep = "_", drop = TRUE)
  
  mm_out <- mm(
    dat,
    choice_clean2 ~ profile1 + Gender + Age + Feminism + Immigration + Crime,
    id = ~ numericalid
  )
  
  mm_out <- mm_out %>%
    mutate(feature = ifelse(feature == "profile1", "Ideology*Immigration", feature)) %>%
    filter(feature == "Ideology*Immigration") %>%
    mutate(
      level = as.character(level),
      pretty = dplyr::case_when(
        level == "Left_No Restrictions"   ~ "Left and No New Restrictions",
        level == "Left_New Restrictions"  ~ "Left and New Restrictions",
        level == "Right_No Restrictions"  ~ "Right and No Restrictions",
        level == "Right_New Restrictions" ~ "Right and New Restrictions",
        TRUE ~ level
      )
    ) %>%
    rename(level_raw = level) %>%
    mutate(level = pretty) %>%
    select(-pretty)
  
  order4 <- c("Left and No New Restrictions",
              "Left and New Restrictions",
              "Right and No Restrictions",
              "Right and New Restrictions")
  present <- intersect(order4, unique(mm_out$level))
  mm_out$level <- factor(mm_out$level, levels = present)
  
  mm_draw(mm_out, title = title, order_levels = present)
}

# Left subgroup (pro‑immigration by respondent attitude)
left_int <- subset(wave2_conjoint,
                   left == 1 & (a6e %in% c("Ni de acuerdo ni en desacuerdo",
                                           "En desacuerdo","Muy en desacuerdo")))
fig6_left  <- mm_interaction(left_int,  "Left Pro-Immigration")
save_fig(file.path(outdir, "mm_left_mig2.png"),  fig6_left,  width = 5, height = 5.4)

# Right subgroup (anti‑immigration by respondent attitude)
right_int <- subset(wave2_conjoint,
                    right == 1 & (a6e %in% c("De acuerdo","Muy de acuerdo")))
fig6_right <- mm_interaction(right_int, "Right Anti-Immigration")
save_fig(file.path(outdir, "mm_right_mig2.png"), fig6_right, width = 5, height = 5.4)


############################################################
# Figure 7: RDD (main effect by election) + Continuity
############################################################

# Loads: d06, d10, d14, d18 (and fields ideo_bin, margin, cep, comuna, etc.)
load(file.path(base_dir, "09_replication/data/rdd_2025aug19.RData"))

# Helper to fit weighted RDD with triangular kernel via felm
run_rdd <- function(dat, yvar = "ideo_bin") {
  bw <- rdbwselect(dat[[yvar]], dat$margin)$bws[1]
  w  <- kernelwts(dat$margin, 0, bw = bw, kernel = "triangular")
  mod <- felm(as.formula(paste0(yvar, " ~ treatment + margin + treatment*margin | cep + comuna | 0 | comuna")),
              data = dat, weights = w, keepCX = TRUE)
  tidy(mod, conf.int = TRUE, conf.level = 0.95)
}

r06 <- run_rdd(d06, "ideo_bin")
r10 <- run_rdd(d10, "ideo_bin")
r14 <- run_rdd(d14, "ideo_bin")
r18 <- run_rdd(d18, "ideo_bin")

pe      <- c(r06$estimate[1], r10$estimate[1], r14$estimate[1], r18$estimate[1])
se      <- c(r06$std.error[1], r10$std.error[1], r14$std.error[1], r18$std.error[1])
upper95 <- c(r06$conf.high[1], r10$conf.high[1], r14$conf.high[1], r18$conf.high[1])
lower95 <- c(r06$conf.low[1],  r10$conf.low[1],  r14$conf.low[1],  r18$conf.low[1])

years <- factor(c("1996-2006", "1996-2010", "1996-2014", "1996-2018"))
results <- data.frame(years, pe, se, lower95, upper95)

fig7_rdd <- ggplot(results, aes(x = years, y = pe)) +
  geom_pointrange(aes(ymin = lower95, ymax = upper95)) +
  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1)) +
  coord_cartesian(ylim = c(-0.2, 1)) +
  scale_x_discrete(breaks = levels(years)) +
  geom_hline(yintercept = 0) +
  xlab("Years") + ylab("RDD estimates") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none",
        text = element_text(size = 20))

save_fig(file.path(outdir, "rdd_ideo_bin_year.png"), fig7_rdd, width = 10, height = 7)


